### This script generates one thousand simulations of the long and short regressions with
### artificial data based on the strata defined in Agnolin Colantone Stanig (2024).


neat.output.models <- array(NA, c(1000,10))

se.coef <- function(model){
  
summary(model)$coefficients[,2]  
  
}

for(m in 1:1000){

      type.A <- data.frame(rep("A",200),rep("lib",200), c(rbinom(100,1,.2 ),rbinom(100,1,.3 )),
                     c(rep("NS", 100), rep("S",100)))

      type.B <- data.frame(rep("B",400),c(rep("lib",200),rep("auth",200)), 
                     c(rbinom(200,1,.3 ),rbinom(200,1,.4 )), c(rep("NS", 200), rep("S",200)))

      type.C <- data.frame(rep("C",200),c(rep("auth",100),rep("auth",100)), c(rbinom(100,1,.7 ),
                                                                        rbinom(100,1,.8 )),
                     c(rep("NS", 100), rep("S",100)))

      names(type.A) <- names(type.B) <- names(type.C) <- c("type", "culture", "RRP", "shocked")

      all.data <- data.frame(rbind(type.A,type.B,type.C))

      all.data$culturenum <- all.data$culture=="auth"

      types <- c(rep("A", 200), rep("B",400),rep("C",200))

      MonY <- lm(RRP~culturenum, data=all.data)

      TonM <- lm(culturenum~shocked, data=all.data)

      TonY <-  lm(RRP~shocked, data=all.data)

      TonYcM <-  lm(RRP~shocked+culturenum, data=all.data)

      neat.output.models[m,] <- c(coef(TonY)[1:2], 
                                se.coef(TonY)[1:2], 
                                coef(TonYcM)[1:3],
                                se.coef(TonYcM)[1:3])

}

save(neat.output.models, file="neat.output.models.RData")

means.sims <- round(colMeans(neat.output.models),2)

tableized <- data.frame(array(NA,c(6,2) ))

#coefficients

tableized[1,1] <- means.sims[1]
tableized[2,1] <- means.sims[2]
tableized[1,2] <- means.sims[5]
tableized[2,2] <- means.sims[6]
tableized[3,2] <- means.sims[7]

#ses
tableized[4,1] <- means.sims[3]
tableized[5,1] <- means.sims[4]
tableized[4,2] <- means.sims[8]
tableized[5,2] <- means.sims[9]
tableized[6,2] <- means.sims[10]

names.coefs <- c("intercept","shock", "culture",paste("se", c("intercept","shock", "culture")))

tableized <- data.frame(names.coefs, tableized)

write.csv(tableized,file="simulations_tableized.csv")

